\documentclass[12pt]{article}
\usepackage{placeins}
\usepackage{amsmath}
\usepackage{Sweave}
\begin{document}
\setkeys{Gin}{width=\textwidth}
\title{Improving the Performance of the HHT using Ensemble Empirical Mode Decomposition}
\author{Daniel Bowman}
\date{\today}
\maketitle
\tableofcontents

\section{Introduction}


This vignette shows how the ensemble empirical mode decomposition method (EEMD) can improve the empirical mode decomposition (EMD) method and the resulting Hilbert spectrogram.
The EEMD is a noise assisted version of the EMD method.
It was developed to reduce aliasing and mode mixing that can occur during EMD.
For more information, see Wu and Huang (2009) ``Ensemble Empirical Mode Decomposition: A Noise Assisted Data Analysis Method.''

This vignette examines three signals.
The first signal demonstrates how the EEMD can improve the quality of intrinsic mode functions (IMFs), particularly at the beginning and end of the signal.

The second signal has two sinusoids with Gaussian white noise added.
This example shows how the EEMD method can improve the IMF set and the resulting Hilbert spectrogram.

The third signal is a transient event recorded at Deception Island volcano, Antarctica.
I show how the EEMD can extract meaningful information from a signal previously thought to be too short for spectral analysis.
I also demonstrate the higher time frequency resolution of the Hilbert spectrogram as opposed to the Fourier spectrogram.

\textbf{WARNING:  The EEMD code in this vignette takes several hours to run.}

\section{Eliminating Spline Fitting Problems}

Consider the following signal:

\begin{equation}
\label{eqn:nonlinsta}
x (t) = \sin(2\pi t+0.5\sin(\pi t))e^{\frac{-(t-50)^{2})}{200}} + \sin(\frac{\pi}{10} t)
\end{equation}

While the EMD successfully removes the nonlinear signal from the lower frequency sinusoid, the first and second IMF contain unusual oscillations in the first 20 and last 20 seconds of the signal (Fig. \ref{fig:nonlinstaemd}).
These oscillations are not implied in the signal equation and likely represent problems with the spline fitting method in the EMD.

If low amplitude Gaussian white noise is added to the signal, the spline fitting issue is solved (Fig. \ref{fig:nonlinstarand}).
However, this creates a lot more IMFs, and the transient nonlinear signal is now distributed across multiple IMFs rather than residing on a single IMF.
This problem is called \emph{mode mixing} and reflects the fact that random noise populates all time scales equally.
Mode mixing can cause aliasing in the Hilbert transform because it disturbs local symmetry in the affected IMFs.

Each trial of the EEMD method adds a different white noise set to the signal and performs EMD. 
If this is done enough times, the averaged IMF set from all the trials approaches the ``true'' signal .
After EEMD, the spline fitting has been improved.
Although the higher frequency signal remains on two IMFs (Fig. \ref{fig:nonlinstasift}), the two averaged modes are smoother compared to the signal with only one EMD run (Fig. \ref{fig:nonlinstarand}).

\subsection{Code}

\textbf{EMD with no added noise}
\begin{Schunk}
\begin{Sinput}
> library(hht)
> dt=0.01
> tt=seq_len(10000)*dt
> sig=sin(2*pi*tt+0.5*sin(pi*tt))*exp(-(tt-50)^2/100)+sin(pi*tt/10)
> emd.result=Sig2IMF(sig, tt)
> PlotIMFs(emd.result)
\end{Sinput}
\end{Schunk}

\textbf{EMD with added noise}
\begin{Schunk}
\begin{Sinput}
> set.seed(628)
> sig2=sin(2*pi*tt+0.5*sin(pi*tt))*exp(-(tt-50)^2/100)+
+     sin(pi*tt/10)+rnorm(length(sig), 0, 0.01)
> emd.result=Sig2IMF(sig2, tt)
> PlotIMFs(emd.result)
\end{Sinput}
\end{Schunk}


\textbf{Perform EEMD}

\begin{Schunk}
\begin{Sinput}
> trials=100
> nimf=11
> noise.amp=0.01
> trials.dir="example_1"
> EEMD(sig2, tt, noise.amp, trials, nimf, 
+     trials.dir = trials.dir, verbose = FALSE)
> EEMD.result=EEMDCompile(trials.dir, trials, nimf)
\end{Sinput}
\end{Schunk}

\textbf{Calculate and Plot Averaged IMFs}

\begin{Schunk}
\begin{Sinput}
> resift.rule="max.var"
> resift.result=EEMDResift(EEMD.result, 
+ resift.rule = resift.rule)
> time.span=c(0, 100)
> imf.list=c(3:7)
> os=TRUE
> res=FALSE
> PlotIMFs(resift.result)
\end{Sinput}
\end{Schunk}

\subsection{Figures}

\FloatBarrier

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-nonlinstaemd}
\end{center}
\caption{EMD of signal in Equation \ref{eqn:nonlinsta}.
A total of 2 IMFs were returned.}
\label{fig:nonlinstaemd}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-nonlinstarand}
\end{center}
\caption{EMD of signal in Equation \ref{eqn:nonlinsta} with random Gaussian noise added.}
\label{fig:nonlinstarand}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-nonlinstasift}
\end{center}
\caption{Major IMFs from EEMD of signal in Equation \ref{eqn:nonlinsta}}.
\label{fig:nonlinstasift}
\end{figure}

\FloatBarrier

\section{Reducing Mode Mixing}

Even a well behaved stationary time series will have problems if there is random noise in addition to the signal.
Fortunately, the EEMD is able to counteract the effect of a noisy signal by averaging a large number of trials together.
As the number of trials grows, the signal tends to emerge from the averaged IMF set.

This example's time series has the form
\begin{equation}
\label{eqn:modemix}
x (t) = 2 \sin(2\pi t) + 0.5\sin(12\pi t) + \mathcal{N}(0, 0.4)
\end{equation}
where 
\begin{equation*}
\mathcal{N}(0, 0.4)
\end{equation*}
produces random Gaussian noise with a mean of 0 and a standard deviation of 0.4.
In theory, the EMD should produce one IMF showing the $2 \sin(2\pi t)$ component and one IMF showing the $0.5\sin(12\pi t)$ component.
However, the presence of noise creates many more IMFs and also causes severe mode mixing.
The most obvious mode mixing occurs where portions of the high frequency sinusoid switch back and forth between IMF 2 and IMF 3.
IMF 4 and IMF 5 both contain portions of the lower frequency, creating severe distortion in that signal band as well (Fig. \ref{fig:modemix}).

The Hilbert spectrogram of this time series should show strong continuous frequency bands at 1 Hz and 6 Hz.
However, the mode mixing produces a chaotic spectrogram with an intermittent low frequency signal around 0.5 Hz and a scattered high frequency signal between 4 and 8 Hz (Fig. \ref{fig:modemixht}).
After a 100 trial EEMD run, the averaged IMF set looks better than the original EMD (Fig. \ref{fig:stackedEEMDimfs}).

The EEMD Hilbert spectrogram is a marked improvement over the original Hilbert spectrogram (Fig. \ref{fig:stackedEEMDht}).
There are two continuous bands of energy that correspond with the sinusoids in Equation \ref{eqn:modemix}.
While the 6 Hz band is still noisy, it is less scattered than in the original spectrogram in Figure \ref{fig:modemixht}.

\subsection{Code}

\textbf{EMD of noisy signal}
\begin{Schunk}
\begin{Sinput}
> library(hht)
> set.seed(628)
> dt=0.01
> tt=seq_len(2000)*dt
> sig=sin(2*pi*tt) + 2*sin(12*pi*tt)+rnorm(length(tt), 0, 0.4)
> emd.result=Sig2IMF(sig, tt)
> PlotIMFs(emd.result)
\end{Sinput}
\end{Schunk}


\textbf{Hilbert spectrogram of noisy signal}
\begin{Schunk}
\begin{Sinput}
> dfreq = 0.05
> freq.span = c(0, 10)
> hgram=HHRender(emd.result, dt, dfreq, freq.span = freq.span, 
+     verbose = FALSE)
> time.span=c(0, 20)
> freq.span=c(0, 10)
> amp.span=c(0.75, 3)
> HHGramImage(hgram, time.span = time.span, freq.span = freq.span,
+     amp.span = amp.span, pretty = TRUE,
+     img.y.format = "%.0f", img.x.format = "%.0f")
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\textbf{Run EEMD}

\begin{Schunk}
\begin{Sinput}
> trials=100
> nimf=9
> noise.amp=0.4
> trials.dir="example_2"
> EEMD(sig, tt, noise.amp, trials, nimf, 
+     trials.dir = trials.dir, verbose = FALSE)
\end{Sinput}
\begin{Soutput}
[1] "Created trial directory: example_2"
\end{Soutput}
\begin{Sinput}
> EEMD.result=EEMDCompile(trials.dir, trials, nimf)
\end{Sinput}
\end{Schunk}

\textbf{Generate averaged IMF set}

\begin{Schunk}
\begin{Sinput}
> PlotIMFs(EEMD.result, imf.list = c(2, 3, 5),
+     res = FALSE, fit.line = TRUE)
\end{Sinput}
\end{Schunk}

\textbf{Render and plot EEMD Hilbert spectrogram}

\begin{Schunk}
\begin{Sinput}
> freq.span=c(0, 10)
> dfreq=0.05
> hgram=HHRender(EEMD.result, dt, dfreq, 
+     freq.span = freq.span, verbose = FALSE)
> time.span=c(0, 20)
> freq.span=c(0, 10)
> amp.span=c(0.75, 3)
> cluster.span=c(5, 100)
> HHGramImage(hgram, time.span = time.span, freq.span = freq.span, 
+     amp.span = amp.span, cluster.span = cluster.span, pretty = TRUE,
+     img.y.format = "%.0f", img.x.format = "%.0f")
\end{Sinput}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}

\subsection{Figures}

\FloatBarrier

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-modemix}
\end{center}
\caption{EMD of signal in Equation \ref{eqn:modemix}.}
\label{fig:modemix}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{eemd_examples-modemixht}
\end{center}
\caption{Hilbert spectrogram of signal in Equation \ref{eqn:modemix}.}
\label{fig:modemixht}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-stackedEEMDimfs}
\end{center}
\caption{Averaged IMF set generated from EEMD of signal in Equation \ref{eqn:modemix}.
Note that the three IMFs shown match the signal content quite well (the red line shows the sum of the IMFs, the black line is the original signal).}
\label{fig:stackedEEMDimfs}
\end{figure}

\begin{figure}[ht]
\begin{center}
\begin{Schunk}
\begin{Soutput}
Adjusting Time and Frequency limits to nice looking numbers (the "pretty" option is currently set to TRUE)
\end{Soutput}
\end{Schunk}
\includegraphics{eemd_examples-stackedEEMDht}
\end{center}
\caption{Hilbert spectrogram generated from EEMD of signal in Equation \ref{eqn:modemix}.}
\label{fig:stackedEEMDht}
\end{figure}

\FloatBarrier

\section{Transient Seismic Event Analyzed with EEMD}

A network of ocean bottom seismometers was deployed in Port Foster  (the flooded caldera of Deception Island volcano) in 2005.
These instruments recorded numerous transient seismic events.
One of these events has been included as an example data set in this R package.
The following example demonstrates how detailed spectral information can be extracted from this signal using the EEMD method.

The EEMD method recovers three main IMFs: a high amplitude, high frequency component, a high amplitude, low frequency component and a low frequency, low amplitude component (Fig. \ref{fig:pferesift}).
While the presence of a high frequency component is clear from the original signal, the spindle shape of this component is not obvious until it is decomposed into an IMF.
Although this could have been discovered by using a bandpass filter, the signal could have been distorted because of the transient nature of the waveform and the fact that this component has frequency modulation.
Therefore, the EEMD has already recovered potentially valuable information on the different components of this signal.

The EEMD Hilbert spectrogram shows that each of the three components varies in frequency throughout the signal (Fig. \ref{fig:pfeht}).
The high frequency component has energy between 15 and 20 Hz, dropping to around 13 Hz in about a quarter of a second.
The low frequency, high amplitude component starts at around 6 Hz, glides to below 5 Hz over a half second, then rises in frequency again for another second or so before tapering off.
The low frequency, low amplitude component decreases from around 5 Hz to around 2 Hz over less than a second.

The Fourier spectrogram cannot resolve such fine details in this signal (Fig. \ref{fig:pfeft}).
It is difficult to determine if there is frequency gliding in the Fourier spectrogram due to the coarseness of the frequency resolution.
This problem could be alleviated by shortening the time window, but that would mask the contribution of lower frequency elements of the signal.
Also, the Fourier spectrogram shows energy arriving slightly before 6.5 seconds, whereas the seismogram shows the main signal starting just after 6.5 seconds.
However, the Hilbert spectrogram shows spectral energy arriving just after 6.5 seconds, thus providing better agreement with the seismogram.

\subsection{Code}

\textbf{Run EEMD}

\begin{Schunk}
\begin{Sinput}
> library(hht)
> data(PortFosterEvent)
> set.seed(628)
> trials=100
> nimf=8
> noise.amp=6.4e-07
> trials.dir="example_3"
> EEMD(sig, tt, noise.amp, trials, nimf, 
+     trials.dir = trials.dir, verbose = FALSE)
\end{Sinput}
\begin{Soutput}
[1] "Created trial directory: example_3"
\end{Soutput}
\begin{Sinput}
> EEMD.result=EEMDCompile(trials.dir, trials, nimf)
\end{Sinput}
\end{Schunk}

\textbf{Create averaged IMF set}

\begin{Schunk}
\begin{Sinput}
> resift.rule="max.var"
> resift.result=EEMDResift(EEMD.result, resift.rule = resift.rule)
> time.span=c(5, 10)
> imf.list=1:3
> os=TRUE
> res=FALSE
> fit.line=TRUE
> PlotIMFs(resift.result, time.span, imf.list, os, res, fit.line)
\end{Sinput}
\end{Schunk}

\textbf{Render and display Hilbert spectrogram}
\begin{Schunk}
\begin{Sinput}
> freq.span=c(0, 30)
> dt = 0.01
> dfreq=0.05
> hgram=HHRender(EEMD.result, dt, dfreq, 
+     freq.span = freq.span, verbose = FALSE)
> time.span=c(5, 9)
> freq.span=c(0, 25)
> amp.span=c(1e-06, 3e-05)
> cluster.span=c(4, 100)
> HHGramImage(hgram, time.span = time.span, freq.span = freq.span,
+     amp.span = amp.span, cluster.span = cluster.span)
\end{Sinput}
\end{Schunk}

\textbf{Render and display Fourier spectrogram}

\begin{Schunk}
\begin{Sinput}
> dt=mean(diff(tt))
> ft=list()
> ft$nfft=4096
> ft$ns=30
> ft$nov=29
> time.span=c(5, 10)
> freq.span=c(0, 25)
> amp.span=c(0.00001, 0.001)
> FTGramImage(sig, dt, ft, time.span = time.span,
+     freq.span = freq.span, amp.span = amp.span)
\end{Sinput}
\end{Schunk}

\subsection{Figures}

\FloatBarrier

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-pferesift}
\end{center}
\caption{Averaged IMFs from EEMD of transient seismic signal.
The red line shows the sum of displayed IMFs; the black line shows the original signal.}
\label{fig:pferesift}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-pfeht}
\end{center}
\caption{Hilbert spectrogram of EEMD of transient seismic signal.}
\label{fig:pfeht}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics{eemd_examples-pfeft}
\end{center}
\caption{Fourier spectrogram of transient seismic signal.}
\label{fig:pfeft}
\end{figure}


\end{document}
